Modelagem geoestatística da profundidade do solo em área de produção silvicultural
Horst
1 de novembro de 2017
1 Objetivo
Realizar a modelagem geoestatística da profundidade do solum área de produção silvicultural
1.1 Caracterização da área de estudo
Os dados que utilizei para este trabalho são provenientes de um povoamento florestal de 108 ha de Pinus taeda L. A área pertence ao município de Campo Belo do Sul, região serrana do Estado de Santa Catarina, Brasil. O clima é do tipo Cfb, mesotérmico, subtropical úmido e com precipitação média de 1.647 mm com chuvas bem distribuídas no ano. A geologia regional é constituída por uma sequência vulcânica de rochas ácidas da Formação Serra Geral, com predomínio de riodacito.
A área possui Neossolos Litólicos e Neossolos Regolíticos, Cambissolos Háplicos e Cambissolos Húmicos, Latossolos Vermelhos e Gleissolos Melânicos, os quais representam a topossequência clássica da área de estudo (Figura 1c). Conforme as tendências observadas no campo, em locais de relevo plano ou suavemente ondulado com boa drenagem estão solos profundos com sequência de horizontes A-Bw (Latossolos), em condições de má drenagem, ocorrem solos com a sequência de horizonte A-Cg (Gleissolos). Nos relevos ondulado ou fortemente ondulado, predominam solos rasos com a sequência de horizontes A ou A-Bi (Neossolos e Cambissolos).
Figure 1.1: Localização da área de estudo no município de Campo Belo do Sul, Estado de Santa Catarina (SC), Brasil (a) e área ampliada com modelo digital de elevação e distribuição dos pontos amostrais (b). Relação da classe de solo e altura das árvores de uma Topossequência típica da área de estudo (c).
1.2 Amostragem
A coleta de dados foi realizada em 94 pontos (Figura 1 b) alocados pelo algoritmo Hipercubo Latino Condicionado (Conditioned Latin Hipercube Sampling – cLHS) através da função clhs do pacote clhs para R. Foram consideradas como variáveis condicionantes da amostragem a elevação, profundidade do vale, índice de umidade topográfico, nivel base da rede de drenagem e declividade as quais, juntas, explicaram aproximadamente 86% da variância topográfica, identificada através da Análise de Componentes Principais (ACP).
Em cada ponto amostral foi medida a profundidade do solum (PF) e a altura total das árvores (h). Consideramos solum a espessura máxima do solo onde as raízes podem se desenvolver sem impedimentos físicos para penetração livre das raízes. Os fatores limitantes considerados foram o lençol freático elevado e o contato com rocha consolidada (contato lítico) com ou sem fissuras.
Para representar a altura das árvores utilizamos o valor médio dos quatro indivíduos de Pinus taeda mais próximos de cada ponto amostral. Portanto, considerando o espaçamento existente entre os indíviuos, cada ponto de amostragem representa uma área de aproximadamente 100 m2 no campo. Apesar do suporte amostral ser areal, considerei-os como sendo em ponto.
Os dados foram armazenados no objeto pontos definido como SpatialPolygonsDataFrame com as coordenadas projetadas em WGS84.
pontos <- read.csv('../data/GateadosDados.csv', dec = ".", sep= ";", stringsAsFactors = FALSE)
sp::coordinates(pontos) <- c('X' , 'Y')
wgs84utm22s <- sp::CRS('+proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
sp::proj4string(pontos) <- wgs84utm22s
pontos <- sp::spTransform(pontos, wgs84utm22s)
# ASR: Não estou entendendo a operação feita aqui. Na linha 73 é definido o WGS 84 zona 22 S como SRC. Em
# seguida, na linha 74 o SRC é transaformado para exatamente o mesmo SCR. Favor verificar.Como uma análise preliminar avaliei a distribuição de frequências dos dados de profundidade do solum. Apesar das amostras apresentarem uma distribuição bimodal (Figura 2), visando atender aos pressupostos teóricos da geoestatística, assumi que os dados provém de uma distribuição normal.
Figure 1.2: Histograma de frequências da profundidade do solum
# Além da distribuição de frequência (acima), é preciso verificar onde estão os dados no espaço.
mapview(pontos, zcol = "PFd")2 Modelagem geoestatística da profundidade do Solum
2.1 Modelo linear misto de variação espacial
Para modelagem geoestatística dos dados foi utilizado o modelo linear misto de variação espacial denotado por
\[Y(\boldsymbol{s}_i) = Z(\boldsymbol{s}_i) + \varepsilon(\boldsymbol{s}_i) = \boldsymbol{x}(\boldsymbol{s}_i)^\text{T}\boldsymbol{\beta} + B(\boldsymbol{s}_i) + \varepsilon(\boldsymbol{s}_i)\]
Para utilizar esse modelo foi necessário supor que os dados são uma realização de um campo aleatório \(Y(\boldsymbol{s}_i)\) com distribuição normal que podem ser descritos como a combinação aditiva de efeitos fixos, efeitos aleatórios e erro aleatório independente.
\(Z(\boldsymbol{s}_i)\) ou sinal possui dois componentes. O primeiro (efeito fixo) \(\boldsymbol{x}(\boldsymbol{s}_i)^\text{T}\boldsymbol{\beta}\) representa os efeitos de origem desterminística, que relaciona a dependência entre a variável e as covariáveis.
O segundo componente do sinal (efeito aleatório), \(B(\boldsymbol{s}_i)\), um campo aleatório Gaussiano estacionário não-observável, descrito por sua função de média e função de covariância.
\(\varepsilon(\boldsymbol{s}_i)\) é o erro (ou ruído), descrito por uma distribuição Gaussina de probabilidade, cujo parâmetro desconhecido de escala é \(\tau\).
2.1.1 Efeito fixo
Defini o efeito fixo do modelo (\(\boldsymbol{x}(\boldsymbol{s}_i)^\text{T}\boldsymbol{\beta}\)) como sendo a profundidade do solum predita via floresta aleatória.
2.1.1.1 Predição via floresta aleatória
As covariáveis utilizadas na predição via floresta aleatória foram covariáveis topográficas e uma covariável de produção.
2.1.1.1.1 Covariáveis topográficas
As covariáveis topográfricas foram derivadas do modelo digital de elevação disponibilizado pelo Governo do Estado de SC - Secretaria de Estado do Desenvolvimento Econômico Sustentável, proveniente do Levantamento Aerofotogramétrico em 2010. Os dados, disponibilizados com resolução de 1 metro, foram reamostrados para resolução espacial de 10 metros utilizando a ferramenta reamostragem no software SAGA GIS.
Os planos de informação utilizados foram importados pela ferramenta raster::rastere armazenados em objetos no formato RasterLayer (Figura 3).
DECLI <- raster::raster("../data/Covars/DECLI.tif")
ELEV <- raster::raster("../data/Covars/ELEV.tif")
TWI <- raster::raster("../data/Covars/TWI.tif")
VD <- raster::raster("../data/Covars/VD.tif")
Figure 2.1: Covariáveis topográficas utilizadas como preditoras no modelo floresta aleatória
A partir da função raster::extract extraí os valores de cada objeto RasterLayer na localização de cada observação contida no objeto espacial pontos.
pontos$DECLI <- raster::extract(DECLI, pontos)
pontos$ELEV <- raster::extract(ELEV, pontos)
pontos$TWI <- raster::extract(TWI, pontos)
pontos$VD <- raster::extract(VD, pontos)2.1.1.1.2 Covariável de produção
Além dos atributos de terreno, foi utilizado um índice de produtividade disponibilizado pela empresa.
Esses índices são utilizados para ordenamento da produção e provém da amostragem de parcelas fixas. A área possui 12 parcelas fixas de inventário contínuo (PIC) de 500 metros quadrados cada que são utilizados na estimativa da produtividade local. Cada PIC é classificada em função da média de altura das 100 árvores de maior perímetro basal da parcela, expressa como altura dominante (Hdom). Em função da Hdom e seus incrementos anuais é estabelecido o índices de sítio (Sitio). Esse índice então é então atribuido à todo o talhão.
As parcelas da área são classificadas em 4 níveis, em que 1 corresponde a melhor e 4 a pior produtividade (Figura 4).
Utilizei a função raster::shapefile para carregar o polígono com informações das sítio - armazenado no objeto ProdutividadePistola. A função sp::spTransform foi usada para projetar as coordenadas original no plano cartesiano (UTM) e a função raster::extract para extraír os valores de cada objeto raster nas localizações de cada observação contidas no objeto espacial pontos.
Sitio <- raster::raster("../data/Covars/Sitio.tif")
Sitio <- as.factor(Sitio)
sp::proj4string(Sitio) <- wgs84utm22s
pontos$Sitio <- raster::extract(Sitio, pontos) %>% as.factor()ProdutividadePistola <-
raster::shapefile('../data/Produtividade/ProdutividadePistola.shp') %>%
sp::spTransform(wgs84utm22s)
ProdutividadePistola$Sitio <- as.factor(ProdutividadePistola$Sitio)
sp::spplot(
ProdutividadePistola, scales = list(draw = T),
main = 'Classificação dos sítios e localização dos pontos') +
lattice::xyplot(Y ~ X, data = as.data.frame(pontos@coords),
pch = 20, col = 'red', lwd = 2, cex = 1.5) %>%
latticeExtra::as.layer()
Figure 2.2: Classificação dos sítios e localização dos pontos na área de estudo
2.1.1.1.3 Construção do modelo de predição e predição espacial
Foi construído um modelo de predição através da função caret::train(method = "rf") que foi armazenado no objeto rf_fit. As variáveis topograficas e a variável de produtividade foram utilizadas como covariávies preditoras. Todas as informações foram utilizadas foram extraídas do objeto espacial pontos.
# Correlação entre covariáveis e profundidade do solo
# Parece não haver correlação alguma (ou muito pequena). Não há correlação com o sítio.
lm(PFd ~ DECLI + ELEV + TWI + VD + Sitio, data = pontos@data) %>% summary()##
## Call:
## lm(formula = PFd ~ DECLI + ELEV + TWI + VD + Sitio, data = pontos@data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.8520 -2.6610 -0.2334 3.0981 5.0357
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.060155 42.297420 0.285 0.7762
## DECLI -16.204984 6.162539 -2.630 0.0101 *
## ELEV -0.002508 0.045630 -0.055 0.9563
## TWI -0.169750 0.339685 -0.500 0.6185
## VD 0.001994 0.052076 0.038 0.9695
## Sitio2 -0.122499 1.133075 -0.108 0.9142
## Sitio3 -1.901585 1.146389 -1.659 0.1008
## Sitio4 -1.304663 1.030864 -1.266 0.2091
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.437 on 86 degrees of freedom
## Multiple R-squared: 0.1493, Adjusted R-squared: 0.08001
## F-statistic: 2.155 on 7 and 86 DF, p-value: 0.04616rf_fit <- caret::train(PFd ~ DECLI + ELEV + Sitio, data = pontos@data, method = "rf", tuneLength = 1, ntree = 1000, importance = TRUE, na.action = na.omit, trControl = trainControl("LOOCV"))
rf_fit## Random Forest
##
## 94 samples
## 3 predictor
##
## No pre-processing
## Resampling: Leave-One-Out Cross-Validation
## Summary of sample sizes: 93, 93, 93, 93, 93, 93, ...
## Resampling results:
##
## RMSE Rsquared
## 3.468303 0.07343996
##
## Tuning parameter 'mtry' was held constant at a value of 2O resultado da predição foi armazenado na coluna rf no objeto espacial pontos e o ajuste entre os dados preditos e observados estão na figura 5.
par <- par(mfrow = c(2,2))
pontos@data$rf <- rf_fit$finalModel$predicted
lm(PFd ~ rf, data = pontos) %>% plot()
Figure 2.3: Ajuste de dados preditos e observados
lm(PFd ~ rf, data = pontos) %>% summary()##
## Call:
## lm(formula = PFd ~ rf, data = pontos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.5646 -3.0334 -0.3804 3.4595 5.7287
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.3135 1.4368 1.610 0.1108
## rf 0.6126 0.2374 2.581 0.0114 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.479 on 92 degrees of freedom
## Multiple R-squared: 0.06752, Adjusted R-squared: 0.05738
## F-statistic: 6.662 on 1 and 92 DF, p-value: 0.01143O modelo rf_fit foi utilizado para fazer a predição espacial em toda área de estudo através da função raster::predict (Figura 6).
beginCluster()
prediction <-
clusterR(brick(DECLI, ELEV, Sitio), raster::predict,
args = list(model = rf_fit, type = "raw", index = 1))
endCluster()
plot(prediction, main = 'Mapa predito de profundidade do solum (dm)')
Figure 2.4: Mapa de profundidade do solum predito pelo modelo floresta aleatória
O RasterLayer Sitio possuia alguns NA (pixel sem valor), logo, a predição não foi realizada nesses pontos. Assim, excluí do RasterLayer prediction todos os pixels sem valor, através da função na.exclude.
O resultado dessa predição espacial foi utilizado como efeito fixo do modelo linear misto de variação.
# Da maneira como está, este comando não possui qualquer efeito.
na.exclude(prediction)## class : RasterLayer
## dimensions : 125, 210, 26250 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 516761.6, 518861.6, 6908574, 6909824 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 2.206938, 8.848382 (min, max)2.2 Variograma amostral
O variograma amostral (Figura 7) foi computado através da função georob::sample.variogram. O estimador para semivariãncia foi Metheron (método dos momentos). Para a obtenção dos parâmetros utilizei um corte de 66% da distância máxima entre os pontos, armazenada no objeto distmax, excluindo os pares de longo alcance.
distmax <-dist(pontos@coords) %>% max() /3
limites <- seq(0, distmax, length.out = 20)
vario<- georob::sample.variogram(PFd ~ rf,
data= pontos, locations = ~ X + Y, lag.dist.def = limites, estimator = "matheron") %>%
plot(ylab = 'Seminvariância', xlab = 'Distância de separação (m)', annotate.npairs = TRUE, main = "Semivariograma")
Figure 2.5: Variograma amostral
2.3 Ajuste da função que descreve o varigrama amostral
O método de ajuste empregado no variograma amostral foi o de quadrados mínimos não-lineares ponderados, com ajuste de um modelo exponencial, com ponderação definida conforme o método de “Cressie”. O processo de estimativa dos parâmetros do modelo exponencial do variograma foi conduzido via otimização usando a função stats::optim(method = "BFGS"). A função resultante desse processo foi armazenada no objeto vario_fit.
vario_fit <-
georob::fit.variogram.model(
vario, variogram.model = 'RMexp', param = c(variance = 10, nugget = 5, scale = 70), weighting.method = "cressie", method = "BFGS")
summary(vario_fit)##
## Call:georob::fit.variogram.model(sv = vario, variogram.model = "RMexp",
## param = c(variance = 10, nugget = 5, scale = 70), weighting.method = "cressie",
## method = "BFGS")
##
## Convergence in 72 function and 29 Jacobian/gradient evaluations
##
## Residual Sum of Squares: 63.64244
##
## Residuals (epsilon):
## Min 1Q Median 3Q Max
## -4.3131 -1.8347 -0.5159 1.1462 3.6763
##
## Variogram: RMexp
## Estimate Lower Upper
## variance 12.69262 9.10587 17.69
## snugget(fixed) 0.00000 NA NA
## nugget 1.33044 0.05434 32.58
## scale 70.37501 50.43150 98.20O ajuste do modelo exponencial ao variograma amostral é mostrado na Figura 8. A curva ajustada passa próximo ao centro de massa dos vinte pontos do variograma amostral.
plot(vario, xlab = 'Distância de separação (m)', ylab = 'Semivariância', annotate.npairs = TRUE)
lines(vario_fit, col = "red", lty = 'dashed')
Figure 2.6: Variograma amostral (em preto) e função exponencial a ele ajustado (vermelho)
Considerando que durante a obteção dos dados de campo as informações de profundidade do solum os valores foram obtidos em decímetros, considerei que a variância do erro de medida é igual a 0,5 dm devido ao arredondamento dos valores.
Assim, foi possível discretizar a variância não explicada em erro de medida corresponde 0,25 do parâmetro nugget. A variância restante foi atribuída a variação espacial não auto-correlacionada espacialmente - não capturada pelo plano amostral snugget. A função resultante desse processo foi armazenado no objeto vario_fit_error.
# Melhor usar um objeto para definir o valor do efeito pepita
nugget <- 0.25
vario_fit_error <- georob::georob(
PFd ~ rf, pontos, locations = ~ X + Y, variogram.model = 'RMexp',
param = c(variance = vario_fit$variogram.object[[1]]$param[['variance']],
nugget = nugget,
snugget = vario_fit$variogram.object[[1]]$param[['nugget']] - nugget,
scale = vario_fit$variogram.object[[1]]$param[['scale']]),
fit.param = georob::default.fit.param(nugget = FALSE, snugget = TRUE),
tuning.psi = 1000, control = georob::control.georob(initial.fixef = 'lm'))
summary(vario_fit_error)##
## Call:georob::georob(formula = PFd ~ rf, data = pontos, locations = ~X +
## Y, variogram.model = "RMexp", param = c(variance = vario_fit$variogram.object[[1]]$param[["variance"]],
## nugget = nugget, snugget = vario_fit$variogram.object[[1]]$param[["nugget"]] -
## nugget, scale = vario_fit$variogram.object[[1]]$param[["scale"]]),
## fit.param = georob::default.fit.param(nugget = FALSE, snugget = TRUE),
## tuning.psi = 1000, control = georob::control.georob(initial.fixef = "lm"))
##
## Tuning constant: 1000
##
## Convergence in 5 function and 3 Jacobian/gradient evaluations
##
## Estimating equations (gradient)
##
## xi scale
## Gradient : -7.950822e-02 -3.915233e-02
##
## Maximized restricted log-likelihood: -246.5921
##
## Predicted latent variable (B):
## Min 1Q Median 3Q Max
## -5.1318 -3.2722 -0.0146 3.9040 4.5422
##
## Residuals (epsilon):
## Min 1Q Median 3Q Max
## -0.257059 -0.064143 -0.002763 0.065310 0.206518
##
## Standardized residuals:
## Min 1Q Median 3Q Max
## -2.38759 -0.81571 -0.03406 0.89906 2.05241
##
##
## Gaussian REML estimates
##
## Variogram: RMexp
## Estimate Lower Upper
## variance 11.518 7.852 16.9
## snugget(fixed) 1.080 NA NA
## nugget(fixed) 0.250 NA NA
## scale 63.738 30.820 131.8
##
##
## Fixed effects coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.8325 1.5463 3.125 0.00238 **
## rf 0.1718 0.2463 0.698 0.48711
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error (sqrt(nugget)): 0.5
##
## Robustness weights:
## All 94 weights are ~= 1.A comparação entre as funções vario_fit_error e vario_fit está na figura 9.
plot(vario)
lines(vario_fit, col = "blue")
lines(vario_fit_error, col = "red")
Figure 2.7: Variograma amostral (em preto) do modelo linear da profundidade do solum e a função exponencial a ele ajustada (azul) e a função ajustada com erro de medida fixo
A função vario_fit_error foi utilizada para a predição espacial. Para isso foi criado um grid de predição grid. A partir da função raster::extract foram extraídos os valores do objeto prediction para o grid.
grid <- sp::spsample(ProdutividadePistola, 10000, type = 'regular')
grid$rf <- raster::extract(prediction, grid)
grid<-
sp::SpatialPointsDataFrame(
coords = grid@coords,
data = data.frame(grid),
proj4string = grid@proj4string)
colnames(grid@coords) <- colnames(pontos@coords)2.4 Predição espacial da resposta
pred_ponto <- predict(
vario_fit_error, newdata = grid, type= "response", signif = 0.95,
control = georob::control.predict.georob(extended.output = TRUE))
sp::gridded(pred_ponto) <- TRUE
spplot(pred_ponto)
Figure 2.8: Mapas de predição - saída extendida do georob
sp::spplot(pred_ponto, zcol = 'pred', main = "Mapa predito de profundidade do solum (dm)")
Figure 2.9: Mapa predito de profundidade do solum (dm)
2.4.1 Cross-validation
validacao <- georob::cv(vario_fit_error, nset = 93)
summary(validacao)##
## Statistics of cross-validation prediction errors
## me mede rmse made qne msse medsse
## -0.02482 -0.11977 3.28538 4.13832 3.42663 1.08653 0.71806
## crps
## 1.911591 - sum((validacao$pred$data - validacao$pred$pred)^2) / sum((validacao$pred$data - mean(validacao$pred$data))^2)## [1] 0.1503402plot(validacao)